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Abstract. The purpose of this review-and-research paper is twofold: (i) to review the role played 
in climate dynamics by fluid-dynamical models; and (ii) to contribute to the understanding and 
reduction of the uncertainties in future climate-change projections. To illustrate the first point, 
we review recent theoretical advances in studying the wind-driven circulation of the oceans. In 
doing so, we concentrate on the large-scale, wind-driven flow of the mid-latitude oceans, which is 
dominated by the presence of a larger, anticyclonic and a smaller, cyclonic gyre. The two gyres 
share the eastward extension of western boundary currents, such as the Gulf Stream or Kuroshio, 
and are induced by the shear in the winds that cross the respective ocean basins. The boundary 
currents and eastward jets carry substantial amounts of heat and momentum, and thus contribute 
in a crucial way to Earth's climate, and to changes therein. 

Changes in this double- gyre circulation occur from year to year and decade to decade. We study 
this low-frequency variability of the wind-driven, double-gyre circulation in mid-latitude ocean 
basins, via the bifurcation sequence that leads from steady states through periodic solutions and 
on to the chaotic, irregular flows documented in the observations. This sequence involves local, 
pitchfork and Hopf bifurcations, as well as global, homoclinic ones. 

The natural climate variability induced by the low-frequency variability of the ocean circulation 
is but one of the causes of uncertainties in climate projections. The range of these uncertainties 
has barely decreased, or even increased, over the last three decades. Another major cause of such 
uncertainties could reside in the structural instability — in the classical, topological sense — of the 
equations governing climate dynamics, including but not restricted to those of atmospheric and 
ocean dynamics. 

We propose a novel approach to understand, and possibly reduce, these uncertainties, based on 
the concepts and methods of random dynamical systems theory. The idea is to compare the climate 
simulations of distinct general circulation models (GCMs) used in climate projections, by applying 
stochastic-conjugacy methods and thus perform a stochastic classification of GCM families. This 
approach is particularly appropriate given recent interest in stochastic parametrization of subgrid- 
scale processes in GCMs. 

As a very first step in this direction, we study the behavior of the Arnol'd family of circle maps 
in the presence of noise. The maps' fine-grained resonant landscape is smoothed by the noise, thus 
permitting their coarse-grained classification. 



1991 Mathematics Subject Classification. 37B55, 37G15, 34F05, 34C37, 34D30, 37C20, 37C15, 37C29, 37E10, 
37E45, 34D08, 37D20, 76U05, 35Q35, 86A05, 86-06. 

Key words and phrases. Climate change, physical oceanography, generic properties of dynamical systems and struc- 
tural stability, random dynamical systems, equations and systems with randomness, bifurcation problems, stochastic 
structural stability, Arnol'd tongues, maps of the circle. 

The present manuscript has been published as: M. Ghil, M.D. Chekroun and E. Simonnet, Physica D: Non- 
linear phenomena. Special Issue on The Euler Equations: 250 Years On, 237(2008), 2111-2126. 

1 



2 



M. GHIL, M.D. CHEKROUN, AND E. SIMONNET 



1. Introduction 

Charney et al. [1] were the first to attempt a consensus estimate of the equilibrium sensitivity 
of climate to changes in atmospheric CO2 concentrations. The result was the now famous range 
for an increase of 1.5 K to 4.5 K in global near-surface air temperatures, given a doubling of CO2 
concentration. 

As the relatively new science of climate dynamics evolved through the 1980s and 1990s, it 
became quite clear — from observational data, both instrumental and paleoclimatic, as well as 
model studies — that Earth's climate never was and is unlikely to ever be in equilibrium. The 
three successive IPCC reports (1991 [2J, 1996, and 2001 [3]) concentrated therefore, in addition to 
estimates of equilibrium sensitivity, on estimates of climate change over the 21st century, based 
on several scenarios of CO2 increase over this time interval, and using up to 18 general circulation 
models (GCMs) in the fourth IPCC Assessment Report (AR4) @]. 

The GCM results of temperature increase over the coming 100 years have stubbornly resisted 
any narrowing of the range of estimates, with results for T s in 2100 as low as 1.4 K or as high 
as 5.8 K, according to the Third Assessment Report. The hope in the research leading up to the 
AR4 was that a set of suitably defined "better GCMs" would exhibit a narrower range of year-2100 
estimates, but this does not seem to have been the case. 

The difficulty in narrowing the range of estimates for either equilibrium sensitivity of climate or 
for end-of-the-century temperatures is clearly connected to the complexity of the climate system, 
the multiplicity and nonlinearity of the processes and feedbacks it contains, and the obstacles to a 
faithful representation of these processes and feedbacks in GCMs. The practice of the science and 
engineering of GCMs over several decades has amply demonstrated that any addition or change in 
the model's "parametrizations" — i.e., of the representation of subgrid-scale processes in terms of 
the model's explicit, large-scale variables — may result in noticeable changes in the model solutions' 
behavior. 

As an illustration, Fig. [T] shows the sensitivity of an atmospheric GCM, which does not include a 
dynamical ocean, to changes in its model parameters. Several thousand simulations were performed 
as part of the "climatepredzctzon.net" experiment [6], using perturbations in several parameters of 
the Hadley Centre's HadAM3 model [7J, coupled to a passive, mixed-layer ocean model. The lower 
panel of Fig. [l] clearly illustrates a wide range of responses to CO2 doubling, from about — 1 K to 
about 8 K |S]. 

The last IPCC report [I] has investigated climate change as a result of various scenarios of CO2 
increase for a set of 18 distinct GCMs. The best estimate of the temperature increase at the end 
of the 21st century from AR4 is about 4.0° C for the worst scenario of greenhouse-gas increase, 
namely A1F1, this scenario envisages, roughly speaking, a future world with a very rapid economic 
growth. The likely range of end-of-century increase in global temperatures is of 2.4-6.4° C in this 
case, and comparably large ranges of uncertainties obtain for all the other scenarios as well [4J. The 
consequences of these scentific uncertainties for the ethical quandaries arising in the socio-economic 
and political decision-making process involved in adaptation to and mitigation of climate changes 
are discussed in [5]. 

An essential contributor to this range of uncertainty is natural climate variability [9] of the 
coupled ocean-atmosphere system. As mentioned already in [10J, most GCM simulations do not 
exhibit the observed interdecadal variability of the oceans' buoyancy- driven, thermohaline circula- 
tion [llj. This circulation corresponds to a slow, pole-to-pole motion of the oceans' main water 
masses, also referred to as the overturning circulation. Cold and denser waters sink in the subpolar 
North Atlantic and lighter waters rise over much wider areas of the lower and southern latitudes. 

Another striking example of low-frequency, interannual-and-interdecadal variability is provided 
by the near-surface, wind-driven ocean circulation [ITJ [12]. Key features of this circulation are 
described at length in Section 2. The influence of strong thermal fronts — like the Gulf Stream in 
the North Atlantic or the Kuroshio in the North Pacific — on the mid-latitude atmosphere above 
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FIGURE 1. Frequency distributions of global mean, annual mean, near-surface tem- 
perature (T g ) for (a) 2,017 GCM simulations, and doubled CO2; and for (b) a subset 
of 414 stable simulations, without substantial climate drift (from [8]). 



is severely underestimated. Typical spatial resolutions in the century-scale GCM simulations of 
[3 El [U d |8j are of the order of 100 km at best, whereas resolutions of 20 km and less would be 
needed to really capture the strong mid-latitude ocean-atmosphere coupling just above the oceanic 
fronts [HI E]. 

An important additional source of uncertainty comes from the difficulty to correctly parametrize 
global and regional effects of clouds and their highly complex small-scale physics. This difficulty is 
particularly critical in the tropics, where large-scale features such as the El-Niho/Southern Oscilla- 
tion and the Madden- Julian oscillation are strongly coupled with convective phenomena [T5 l [T6 } [T7] . 

The purpose of this paper is twofold. First, we describe in Section 2 the most recent theoretical 
results regarding the internal variability of the mid-latitude wind-driven circulation, viewed as a 
problem in nonlinear fluid mechanics. These results rely to a large extent on the deterministic 
theory of dynamical systems [HI [19]. Second, we address in Section 3 the more general issue of 
uncertainties in climate change projections. Here we rely on concepts and methods from random 
dynamical systems theory [20] to help understand and possibly reduce these uncertainties. Much 
of the material in the latter section is new; it is supplemented by rigorous mathematical definitions 
and results in Appendices A and B. A summary and an outlook on future work follow in Section 4. 



2. Natural variability of the wind-driven ocean circulation 

2.1. Observations. To a first approximation, the main near-surface currents in the oceans are 
driven by the mean effect of the winds. The trade winds near the equator blow mainly from east to 
west and are called also the tropical easterlies. In mid-latitudes, the dominant winds are the pre- 
vailing westerlies, and towards the poles the winds are easterly again. Three of the strongest near- 
surface, mid-and-high-latitude currents are the Antarctic Circumpolar Current, the Gulf Stream 
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in the North Atlantic, and the Kuroshio Extension off Japan. The Antarctic Circumpolar Current, 
sometimes called the Westwind Drift, circles eastward around Antarctica; see Fig. 2. 




Warm Current Cold Current 



Figure 2. A map of the main oceanic currents: warm currents in red and cold ones 
in blue. 

The Gulf Stream is an oceanic jet with a strong influence on the climate of eastern North America 
and of western Europe. Actually, the Gulf Stream is part of a larger, gyre-like current system, which 
includes the North Atlantic Drift, the Canary Current and the North Equatorial Current. It is also 
coupled with the pole-to-pole overturning circulation. From Mexico's Yucatan Peninsula, the Gulf 
Stream flows north through the Florida Straits and along the East Coast of the United States. Near 
Cape Hatteras, it detaches from the coast and begins to drift off into the North Atlantic towards 
the Grand Banks near Newfoundland. 

The Coriolis force is responsible for the so-called Ekman transport, which deflects water masses 
orthogonally to the near-surface wind direction and to the right [2TJ |22l [23] . In the North At- 
lantic, this Ekman transport creates a divergence and a convergence of near-surface water masses, 
respectively, resulting in the formation of two oceanic gyres: a smaller, cyclonic one in subpolar 
latitudes, the other larger and anticyclonic in the subtropics. This type of double-gyre circulation 
characterizes all mid-latitude ocean basins, including the South Atlantic, as well as the North and 
South Pacific. 

The double-gyre circulation is intensified as the currents approach the East Coast of North 
America due to the /3-effect. This effect arises primarily from the variation of the Coriolis force 
with latitude, while the oceans' bottom topography also contributes to it. The former, planetary 
/3-effect is of crucial importance in geophysical flows and induces free Rossby waves propagating 
westward [2T1 [22| [23] . 

The currents along the western shores of the North Atlantic and of the other mid-latitude ocean 
basins exhibit boundary-layer characteristics and are commonly called western boundary currents 
(WBCs). The northward- flowing Gulf Stream and the southward- flowing Labrador Current exten- 
sion meet near Cape Hatteras and yield a strong eastward jet. The formation of this jet and of the 
intense recirculation vortices near the western boundary, to either side of the jet, is mostly driven 
by internal, nonlinear effects. 

Figure [3] illustrates how these large-scale wind-driven oceanic flows self-organize, as well as 
the resulting eastward jet. Different spatial and time scales contribute to this self-organization, 
mesoscales eddies playing the role of the synoptic-scale systems in the atmosphere. Warm and cold 
rings last for several months up to a year and have a size of about 100 km; two cold rings are 
clearly visible in Fig. [3| Meanders involve larger spatial scales, up to 1000 km, and are associated 
with interannual variability. The characteristic scale of the jet and gyres is of several thousand 
kilometers and they exhibit their own intrinsic dynamics on time scales of several years to possibly 
several decades. 

A striking feature of the wind-driven circulation is the existence of two well-known North- Atlantic 
oscillations, with a period of about 7 and 14 years, respectively. Data analysis of various climatic 
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Figure 3. A satellite image of the sea surface temperature (SST) over the north- 
western North Atlantic (U.S. National Oceanic and Atmospheric Administration), 
together with a sketch of the associated double-gyre circulation. An idealized view 
of the amount of potential vorticity injected into the ocean circulation by the trade 
winds, westerlies and polar easterlies is shown to the right. 



variables, such as sea surface temperature (SST) over the North Atlantic or sea level pressure (SLP) 
over western Europe [2U|25l|26] and local surface air temperatures in Central England [27], as well 
as of proxy records, such as tree rings in Britain, travertine concretions in southeastern France [28], 
and Nile floods over the last millennium or so [29] , all exhibit strikingly robust oscillatory behavior 
with a 7-yr period and, to a lesser extent, with a 14-yr period. Variations in the path and intensity 
of the Gulf Stream are most likely to exert a major influence on the climate in this part of the world 
[30] . This is why theoretical studies of the low- frequency variability of the double-gyre circulation 
are important. 

Given the complexity of the processes involved, climate studies have been most successful when 
using not just a single model but a full hierarchy of models, from the simplest "toy" models to the 
most detailed GCMs [T7] . In the following, we describe one of the simplest models of the hierarchy 
used in studying this problem. 



2.2. A simple model of the double- gyre circulation. The simplest model that includes many 
of the mechanisms described above is governed by the barotropic quasi- geostrophic (QG) equations. 
The term geostrophic refers to the fact that large-scale rotating flows tend to run parallel to, rather 
than perpendicular to constant-pressure contours; in the oceans, these contours are associated with 
the deviation from rest of the surfaces of equal water mass, due to Ekman pumping. Geostrophic 
balance implies in particular that the flow is divergence-free. The term barotropic, as opposed to 
baroclinic, has a slightly different meaning in geophysical fluid dynamics than in engineering fluid 
mechanics: it means that the model describes a single fluid layer of constant density and therefore 
the solutions do not depend on depth [211 [23] . 

We consider an idealized, rectangular basin geometry and simplified forcing that mimics the 
distribution of vorticity contribution by the winds, as sketched to the right of Fig. [3| In our 
idealized model, the amounts of subpolar and subtropical vorticity injected into the basin are equal 
and the rectangular domain Q = (0,L X ) x (0, L y ) is symmetric about the axis of zero wind stress 
curl. The barotropic two-dimensional (2-D) QG equations in this idealized setting are: 

q t + J 00, q) - vA 2 ^ + /xA^ = -r sin ^ , 

(2.1) 

q = Aip - XriP + /3y. 
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Here q and i/j are the potential vorticity and streamfunction, respectively, and the Jacobian J 
corresponds to the advection of potential vorticity by the flow, J(^,g) = ifj x Qy — ^ y q x — u • Vg, 
where u = (— t/^,^), x points east and y points north. The physical parameters are the strength 
of the planetary vorticity gradient /3, the Rossby radius of deformation A^ 2 , the eddy- viscosity 
coefficient is, the bottom friction coefficient fi, and the wind-stress intensity r. We use here free- 
slip boundary conditions ij) — A 2 ^ = 0; the qualitative results described below do not depend on 
the particular choice of homogeneous boundary conditions. 



We consider (2.1) as an infinite-dimensional dynamical system and study its bifurcation sets as 



the parameters change. Two key parameters are the wind stress intensity r and the eddy viscosity 



v. An important property of (2.1) is its mirror symmetry in the y — L y /2 axis. This symmetry 
can be expressed as invariance with respect to the discrete Z2 group S: 

(2.2) S[^(x,y)] = -^(x,L y -y); 



any solution of (2.1) is thus accompanied by its mirror-conjugated solution. Hence, in generic 
terms, the prevailing bifurcations are of either the symmetry-breaking or the saddle-node or the 
Hopf type. 

2.3. Bifurcations in the double-gyre problem. The historical development of a comprehensive 
nonlinear theory of the double-gyre circulation is interesting on its own, having seen substantial 
progress in the last 15 years. One can distinguish four main steps. 

2.3.1. Symmetry-breaking bifurcations. The first step was to realize that the first generic bifurcation 
of this QG model was a genuine pitchfork bifurcation that breaks the system's symmetry as the 
nonlinearity becomes large enough [3H E21 [33]. The situation is shown in Fig. |] When the 
forcing is weak or the dissipation is large, there is only one steady solution, which is antisymmetric 
with respect to the mid-axis of the basin. This solution exhibits two large gyres, along with 
their typical, /3-induced WBCs. Away from the western boundary, such a near-linear solution 
(not shown) is dominated by Sverdrup balance between wind stress curl and the meridional mass 
transport pTll33] . 

As the wind stress increases, the near-linear Sverdrup solution develops an eastward jet along the 
mid-axis, which penetrates farther into the domain. This more intense, and hence more nonlinear 
solution is still antisymmetric about the mid-axis, but loses its stability for some critical value of 
the wind-stress intensity (indicated by "Pitchfork" in Fig. [4]). 

A pair of mirror-symmetric solutions emerges and is characterized by a rather different vorticity 
distribution; the streamfunction fields associated with the two stable steady-state branches are 
plotted to the upper-left and right of Fig. |4| In particular, the jet in such a solution exhibits a 
large meander, reminiscent of the one seen in Fig. [3]just downstream of Cape Hatteras; note that 
the colors in Fig. [4] have been chosen to facilitate the comparison with Fig. [3j These asymmetric 
flows are characterized by one gyre being stronger in intensity than the other and therefore the jet 
is deflected either to the southeast or to the northeast. 

2.3.2. Gyre modes. The next step was taken in part concurrently with [3TJ [32] and in part shortly 
after (35j [36J |37] the first one. It involved the study of time-periodic instabilities through Hopf 
bifurcation from either an antisymmetric or an asymmetric steady flow. Some of these studies 
concentrated on the wind-driven circulation formulated for the stand-alone, single gyre [371 EH] . 
The idea was to develop a full generic picture of the time-dependent behavior of the solutions in 
more turbulent regimes, by classifying the various instabilities in a comprehensive way. However, it 
quickly appeared that one kind of asymmetric instabilities, called gyre modes [32l[35], was prevalent 
across the full hierarchy of models of the double-gyre circulation; furthermore, these instabilities 
trigger the lowest nonzero frequency present in these models. 

These modes always appear after the first pitchfork bifurcation, and it took several years to really 
understand their genesis: gyre modes arise as two eigenvalues merge — one is associated with a 
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Figure 4. Generic bifurcation diagram for the barotropic QG model of the double- 
gyre problem: the asymmetry of the solution is plotted versus the intensity of the 
wind stress r. The streamfunction field is plotted for a steady-state solution asso- 
ciated with each of the three branches; positive values in red and negative ones in 
blue (after 



symmetric eigenfunction and responsible for the pitchfork bifurcation, the other is associated with 
an antisymmetric eigenfunction [39]; this merging is marked by M in Fig. [4j 

Such a phenomenon is not a bifurcation stricto sensu: one has topological C° equivalence before 
and after the eigenvalue merging, but not from the C 1 point of view. We recall here that functions 
are C k if they and their inverses are k times continuously differentiate. Still, this phenomenon 
is quite common in small-dimensional dynamical systems with symmetry, as exemplified by the 
unfolding of codimension-2 bifurcations of Bogdanov-Takens type [19J. In particular, the fact that 
gyre modes trigger the lowest- frequency of the model is due to the frequency of these modes growing 
quadratically from zero until nonlinear saturation. Of course, these modes, in turn, become unstable 
shortly after the merging, through a Hopf bifurcation, indicated by "Hopf ' in Fig. [4j 

2.3.3. Global bifurcations. The importance of these gyre modes was further confirmed recently 
through an even more puzzling discovery. Several authors realized, independently of each other, that 
the low-frequency dynamics of their respective double-gyre models was driven by intense relaxation 
oscillations of the jet [IUJ EJ H21 S3 HU SSI SB]- These relaxation oscillations, already described 
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in [32l|35], were now attributed to homoclinic bifurcations, with a global character in phase space 
jT9l [22] . In effect, the QG model reviewed here undergoes a genuine homoclinic bifurcation (see Fig. 
[4]), which is generic across the full hierarchy of double-gyre models. Moreover, this global bifurcation 
is associated with chaotic behavior of the flow due to the Shilnikov phenomenon [JHl 06] , which 
induces horseshoes in phase space. 

The connection between such homoclinic bifurcations and gyre modes was not immediately 
obvious, but Simonnet et al [16] emphasized that the two were part of a single, global dynamical 
phenomenon. The homoclinic bifurcation indeed results from the unfolding of the gyre modes' limit 
cycles. This familiar dynamical scenario is again well illustrated by the unfolding of a codimension- 
2 Bogdanov-Takens bifurcation, where the homoclinic orbits emerge naturally. We deal, once 
more, with the lowest-frequency modes, since homoclinic orbits have an infinite period. Due to 
the genericity of this phenomenon, it was natural to hypothesize that the gyre-mode mechanism, 
in this broader, global-bifurcation context, gave rise to the observed 7-yr and 14-yr North- Atlantic 
oscillations. Although this hypothesis may appear a little farfetched, in view of the simplicity of 
the double-gyre models analyzed in detail so far, it poses an interesting question. 

2.3.4. Quantization and open questions. The chaotic dynamics observed in the QG models after the 
homoclinic bifurcation is eventually destroyed as the nonlinearity and the resolution both increase. 
As one expects the real oceans to be in a far more turbulent regime than those studied so far, some 
authors proposed different mechanisms for low-frequency variability in fully turbulent flow regimes 
[13 HE]. It turns out, though, that — just as gyre modes could be reconciled with homoclinic- 
driven dynamics, — the latter can also be reconciled with eddy-driven dynamics, via the so-called 
quantization of the low- frequency dynamics [^9] . 

Primeau [50] showed that, in large basins comparable in size with the North Atlantic, there is 
not only one but a set of successive pitchfork bifurcations. One supercritical pitchfork bifurcation, 
associated with the destabilization of antisymmetric flows, is followed generically by a subcritical 
one, associated this time with a stabilization of antisymmetric flows (modulo high-frequency insta- 
bilities) [49J. As a matter of fact, this phenomenon appears to be a consequence of the spectral 
behavior of the 2-D Euler equations [51], and hence of the closely related barotropic QG model in 
bounded domains. 

Remarkably, this scenario repeats itself as the nonlinearity increases, but now higher wavenum- 
bers are involved in physical space. Simonnet [49] showed that this was also the case for gyre 
modes and the corresponding dynamics induced by global bifurcations: the low-frequency dynam- 
ics is quantized as the jet stream extends further eastward into the basin, due to the increased 
forcing and nonlinearity. Figure [5] illustrates this situation: two families of regimes can be identi- 
fied, the colored bands correspond to (supercritical) regimes driven by the gyre modes, the others 
to (subcritical) regimes driven by the eddies. Note that this scenario is also robust to perturbing 
the problem's symmetry. 

The successive-bifurcation theory appears therewith to be fairly complete for barotropic, single- 
layer models of the double-gyre circulation. This theory also provides a self-consistent, plausible 
explanation for the climatically important 7-year and 14-year oscillations of the oceanic circulation 
and the related atmospheric phenomena in and around the North- Atlantic basin [HI [12j [2H |25l [26J, 
E2 EHl [291 HSl HE] . The dominant 7- and 14-year modes of this theory also survive perturbation by 
seasonal-cycle changes in the intensity and meridional position of the westerly winds [52J. 

In baroclinic models, with two or more active layers of different density, baroclinic instabilities 
[HI E31 12H 1221 1231 IBQ1 13HI SSI S3 SS] surely play a fundamental role, as they do in the observed 
dynamics of the oceans. However, it is not known to what extent baroclinic instabilities can destroy 
gyre-mode dynamics. The difficulty lies in a deeper understanding of the so-called rectification 
process [53], which arises from the nonzero mean effect of the baroclinic component of the flow. 

Roughly speaking, rectification drives the dynamics far away from any steady states. In this 
situation, dynamical systems theory cannot be used as an explanation of complex, observed behavior 
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Figure 5. Two-parameter plane, with the wind-stress intensity r vs. the eddy- 
viscosity coefficient v\ the curves indicate the locations of supercritical and subcrit- 
ical pitchfork bifurcations. Each band is associated with a different wavenumber 
and timescale (from [49J). 

resulting from successive bifurcations that are rooted in a simple steady state. Other tools from 
statistical mechanics and nonequilibrium thermodynamics should, therefore, be considered [5 ^ (55 } 
l56l l57] . Combining these tools with those of the successive-bifurcation approach may eventually 
lead to a more general and complete physical characterization of gyre modes in realistic models. 

3. CLIMATE-CHANGE PROJECTIONS AND RANDOM DYNAMICAL SYSTEMS (RDSs) 

As discussed in Section 1, the climate system's natural variability and the difficulties in parametriz- 
ing subgrid-scale processes are not the only causes for the uncertainties in projecting future climate 
evolution. In this section, we address more generally these uncertainties and present a novel ap- 
proach for treating them. To do so, we start with some simple ideas about deterministic vs. 
stochastic modeling. 

3.1. Background and motivation. Many physical phenomena can be modeled by deterministic 
evolution equations. Dynamical systems theory is essentially a geometric approach for studying 
the asymptotic, long-term properties of solutions to such equations in phase space. Pioneered by 
H. Poincare [58J, this theory took great strides over the last fifty years. To apply the theory in 
a reliable manner to a set of complex physical phenomena, one needs a criterion to evaluate the 
robustness of a given model within a class of dynamical systems. Such a criterion should help us 
deal with the inescapable uncertainties in model formulation, whether due to incomplete knowledge 
of the governing laws or inaccuracies in determining model parameters. 

In this context, Andronov and Pontryagin |59j took a major step toward classifying dynamical 
systems, by introducing the concept of structural stability. Structural stability means that a small, 
continuous perturbation of a given system preserves its dynamics up to a homeomorphism, i.e., up 
to a one-to-one continuous change of variables that transforms the phase portrait of our system 
into that of the nearby system; thus fixed points go into fixed points, limit cycles into limit cycles, 
etc. Closely related is the notion of hyperbolicity introduced by Smale [60J. A system is hyperbolic 
if, (very) loosely speaking, its limit set can be continuously decomposed into invariant sets that are 
either contracting or expanding; see |61| for more rigorous definitions. 

A very simple example is the phase portrait in the neighborhood of a fixed point of saddle type. 
In this case, the Hartman-Grobman theorem states that the dynamics in this neighborhood is 
structurally stable. The converse statement, i.e. whether structural stability implies hyperbolicity, 
is still an open question; the equivalence between structural stability and hyperbolicity has only 
been shown in the C 1 case, under certain technical conditions [621 [631 EH ES] . Bifurcation theory is 
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well grounded in the setting of hyperbolic dynamics. Problems with hyperbolicity and bifurcations 
arise, however, when one deals with more complicated limit sets. 

Hyperbolicity was introduced initially to help pursue the "dynamicist's dream" of finding, in the 
abstract space of all possible dynamical systems, an open and dense set consisting of structurally 
stable ones. Being open and dense, roughly speaking, means that any possible dynamical system 
can be approximated by systems taken from this set, while systems in its complement are negligible 
in a suitable sense. 

Smale conjectured that hyperbolic systems form an open and dense set in the space of all C 1 
dynamical systems. If this conjecture were true then hyperbolicity would be typical of all dynamics. 
Unfortunately, though, this conjecture is only true for one-dimensional dynamics and flows on disks 
and surfaces [66J. Smale [67] himself found several counterexamples to his conjecture. Newhouse 
|68j was able to generate open sets of nonhyperbolic diffeomorphisms using homoclinic tangencies. 
For the physicist, it is even more striking that the famous Lorenz attractor |69j is structurally 
unstable. Families of Lorenz attractors, classified by topological type, are not even countable 
[70j[7l]. In each of these examples, we observe chaotic behavior in a nonhyperbolic situation, i.e., 
nonhyperbolic chaos. 

Nonhyperbolic chaos appears, therefore, to be a severe obstacle to any "easy" classification of 
dynamic behavior. As mentioned by Palis [65J, Kolmogorov already suggested at the end of the 
sixties that "the global study of dynamical systems could not go very far without the use of new 
additional mathematical tools, like probabilistic ones." Once more, Kolmogorov showed prophetic 
insight, and nowadays the concept of stochastic stability is an important tool in the study of gener- 
icity and robustness for dynamical systems. To replace the failed program of classifying dynamical 
systems based on structural stability and hyperbolicity, Palis [65J formulated the following global 
conjecture: systems having only finitely many attractors (i.e. periodic or chaotic sinks) - such that 
(i) the union of their basins has full Lebesgue measure; and (ii) each is stochastically stable in their 
basins of attraction - are dense in the C r ,r > 1 topology. A system is stochastically stable if its 
Sinai-Ruelle-Bowen (SRB) measure [72J is stable with respect to stochastic perturbations, and the 
SRB measure is given by lim n ^ 00 ^ ^2 i S Zi , with zi being the successive iterates of the dynamics. 
This measure is obtained intuitively by allowing the entire phase space to flow onto the attractor 

Stochastic stability is fundamentally based on ergodic theory. We would like to consider a more 
geometric approach, which can provide a coarser, more robust classification of GCMs and their 
climate-change projections. In this section, we propose such an approach, based on concepts from 
the rapidly growing field of random dynamical systems (RDSs), as developed by L. Arnold [20] 
and his "Bremen group," among others. RDS theory describes the behavior of dynamical systems 
subject to external stochastic forcing; its tools have been developed to help study the geometric 
properties of stochastic differential equations (SDEs). In some sense, RDS theory is the stochastic 
counterpart of the geometric theory of ordinary differential equations (ODEs). This approach 
provides a rigorous mathematical framework for a stochastic form of robustness, while the more 
traditional, topological concepts do not seem to be appropriate. 

3.2. RDSs, random attractors, and robust classification. Stochastic parametrizations for 
GCMs aim at compensating for our lack of detailed knowledge on small spatial scales in the best 
way possible [THEE EE HUES El]- The underlying assumption is that the associated time scales are 
also much shorter than the scales of interest and, therefore, the lag correlation of the phenomena 
being parametrized is negligibly small. Stochastic parametrizations thus essentially transform a 
deterministic autonomous system into a nonautonomous one, subject to random forcing. 

Explicit time dependence in a dynamical system immediately raises a technical difficulty. Indeed, 
the classical notion of attractor is not always relevant, since any object in phase space is "moving" 
with time and the natural concept of forward asymptotics is meaningless. One needs therefore 
another notion of attractor. In the deterministic nonautonomous framework, the appropriate notion 
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is that of a pullback attractor [80], which we present below. The closely related notion of random 
attractor in the stochastic framework is also explained briefly below, with further details given in 
Appendix A. 

3.2.1. Framework and objectives. Before defining the notion of pullback attractor, let us recall some 
basic facts about nonautonomous dynamical systems. Consider the ODE 



on a vector space X; this space could even be infinite-dimensional, if we were dealing with partial or 
functional differential equations, as is often the case in fluid-flow and climate problems. Rigorously 
speaking, we cannot associate a dynamical system acting on X with a nonautonomous ODE; 
nevertheless, in the case of unique solvability of the initial- value problem, we can introduce a two- 
parameter family of operators {S(t, s)} t > s acting on X, with s and t real, such that S(t, s)x(s) = 
x(t) for t > «s, where x(t) is the solution of the Cauchy problem with initial data x(s). This family 
of operators satisfies S(s,s) = Idx and S(t,r) o S(r,s) = S(t,s) for all t > r > s, and all real 
s. This family of operators is called a "process" by Sell [81] .It extends the classical notion of the 
resolvent of a nonautonomous linear ODE to the nonlinear setting. 

We can now define the pullback attractor as simply the family of invariant sets {*4(£)} that 
satisfy for every real t and all xq in X: 



"Pullback" attraction does not involve running time backwards; it corresponds instead to the idea 
of measurements being performed at present time t in an experiment that was started at some time 
s < t in the past: the experiment has been running for long enough, and we are thus looking now 
at an "attracting state." Note that there exists several ways of defining a pullback attractor — the 
one retained here is a local one (cf. [80J and references therein); see [82] for further information on 
nonautonomous dynamical systems in general. 

In the stochastic context, noise forcing is modeled by a stationary stochastic process. If the 
deterministic dynamical system of interest is coupled to this stochastic process in a reasonable way 
— to be expressed below by the "cocycle property" — then random pullback attractors may appear. 
These pullback attractors will exist for almost each sample path of the driving stochastic process, so 
that the same probability distribution governs both sample paths and their corresponding pullback 
attractors. A more detailed explanation is given in Appendix A. 

Roughly speaking, this concept of random attractor provides a geometric framework for the 
description of asymptotic regimes in the context of stochastic dynamics. To compare different 
stochastic systems in terms of their random attractors that evolve in time, it would be nice to be 
able to identify the common underlying geometric structures via a random change of variables. 
This identification is achieved through the concept of stochastic equivalence that is developed in 
Appendix A, and it is central in obtaining a coarser and more robust classification than in the 
purely deterministic context. 

Returning now to our main objective, suppose for instance that one is presented with results 
from two distinct GCMs, say two probability distributions functions (PDFs) of the temperature 
or precipitation in a given area. These two PDFs are generated, typically, by an ensemble of each 
GCM's simulations, as described in the introduction, and they are likely to differ in their spatial 
pattern. To ascertain the physical significance of this discrepancy, one needs to know how each 
GCM result varies as either a parametrization or a parameter value are changed. 

In order to consider the difficult question of why GCM responses to CO2 doubling might differ, 
one idea is to investigate the structure of the space of all GCMs. We mean therewith the space 
of all deterministic GCMs, when their stochastic parametrizations are switched off. We know, by 
now, from experience with GCM results over several decades — including the four IPCC assessment 
reports J2j H] and the climateprerfzctzon.net exercise [HI El [8] — that there is enormous scatter in 



(3.1) 



x = f(t,x) 



(3.2) 



lim dist (S(t,s)x ,A(t)) = 0. 
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this space; see also [HUES]. Our question, therefore, is: can we achieve a more robust classification 
of GCMs when stochastic parametrizations are used and for a given level of the noise? 

As mentioned in Section |3.1[ such a classification is not feasible by restricting ourself to de- 
terministic systems and topological concepts. As one switches on stochastic parametrizations 
[7U EH [76l [771 EH1 [79], the situation might change, and hopefully improve, dramatically: as the 
noise level becomes large enough, the models' deterministic behavior may be completely destroyed, 
and all the results could cluster into one huge, diffuse clump. We would like, therefore, to inves- 
tigate how a classification based on stochastic equivalence evolves as the level of the noise or the 
stochastic parametrizations change. As the noise tends to zero, do we recover the "granularity" of 
the set of all deterministic dynamical systems? This idea is schematically represented in Fig. |6j 
for a given level of the noise, we expect the space of all GCMs to be decomposed into a possibly 
finite number of classes. Within one of these classes, all the GCMs are topologically equivalent in 



the stochastic sense defined above; see Eq. (A. 2) 




Figure 6. A conjectural view of stochastic classification for GCMs, using the con- 
cept of random attractors. Each point in red represents a GCM in which stochas- 
tic parametrizations are switched off, while each gray area represents a cluster of 
stochastically equivalent GCMs for a given level of the noise. 

Serious difficulties might arise in this program, due to the presence of nonhyperbolic chaos in 
climate models. Several studies have pointed out that the characteristics of nonhyperbolic chaos 
in the presence of noise may depend on its intensity and statistics [86l [871 188| 189]- 

Such issues, however, go well beyond the setting of this paper and are left for further investigation. 
Much more modestly, we will study here whether, in certain very simple cases, the conjectural view 
of Fig. [6] might be relevant for some dynamical systems that are "metaphors" of climate dynamics. 
The following subsection is dedicated to the study of such a metaphorical object, namely the 
Arnol'd circle map. 

3.2.2. The stochastically perturbed circle map. To go beyond our pictorial view of stochastic classi- 
fication for GCMs in Fig. [6j we study now the effect of noise on a family of diffeomorphisms of the 
circle. This toy model exhibits two features of interest for our purpose. The first one is that the 



two-parameter family {F T ^ e } defined by Eq. (3.3) below exhibits an infinite number of topological 



classes |18| . The second feature of interest is the frequency-locking behavior observed in many 
field of physics in general [90l EU [92] and in some El-Niho/Southern-Oscillation (ENSO) models in 
particular [931 EH [951 ESI E3 EHl EH]. Studying noise effects on these two features has, therefore, 
physical and mathematical, as well as climatological relevance. 

Many physical and biological systems exhibit interference effects due to competing periodici- 
ties. One such effect is mode locking, which is due to nonlinear interaction between an "internal" 
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frequency uj{ of the system and an "external" frequency ou e . In the ENSO case, the external peri- 
odicity is the seasonal cycle. A simple model for systems with two competing periodicities is the 
well-known Arnol'd family of circle maps 



x n+1 = F T ^(x n ) := x n + r - esin(27rx n ) mod 1, 

e and e parameterizes the magnitude of nonlinear effects; the map (|3.3D 



(3.3) 

where basically r := uji/uj, 
is often called the standard circle map [T8] . 

These maps also represent frequency locking near a bifurcation of Neimark-S acker type {e.g. 
[100] . p. 434); here the parameter r is typically interpreted as the novel (internal) frequency 
involved in the bifurcation and e corresponds to the nonlinearity near the bifurcation. 

Such nonlinear coupling between two oscillators gives rise to a characteristic pattern, in the plane 
of e vs. r, called Arnol'd tongues. We computed this pattern numerically for the family of Eq. 
(3.3), together with a cross-section at a fixed value of e; see Fig. [7| This cross-section exhibits the 
so-called Devil's staircase, with "steps" on which the rotation number [58J is constant within each 
Arnol'd tongue; the rotation number measures the average rotation per iterate of (3.3). 




200 400 



Figure 7. 
r and e are 
the right. 



Arnol'd tongues for the family of diffeomorphisms of the circle; units for 
5 • 10 -4 and 10 -4 respectively. Devil's staircase in the cross-section to 



For e = 0, two types of phenomena occur: either r is rational and in this case the dynamics 
is periodic with period g, where r = p/q, or r is irrational and the iterates {x n } fill the whole 
circle densely. As e increases, an Arnol'd tongue of increasing width grows out of each r = p/q on 
the abscissa e = 0. It follows that, in this very simple case, such an Arnol'd tongue corresponds 
to hyperbolic dynamics that is robust to perturbations, as verified by linearizing the map at the 
periodic point; the rotation number is then rational and equal to p/q. 

The set of all these tongues is dense within the whole circle map family, while the Lebesgue 
measure of this set, at given e, tends to zero as e goes to zero. On the contrary, if a point in the 
(t, e)-plane does not belong to an Arnol'd tongue, the rotation number for those parameter values 
is irrational and the dynamics is nonhyperbolic; the latter fact follows, for instance, from a theorem 
of Denjoy |101j showing that such dynamics is smoothly equivalent to an irrational rotation. The 
probability to observe nonhyperbolic dynamics tends therewith to unity as e goes to zero. One has, 
therefore, a countably infinite number of distinct topological classes, namely the Arnol'd tongues 
p/q, and an uncoutably infinite number of maps with irrational rotation numbers. 

What happens when noise is added in Eq. (3.3)? We consider here the case of additive forcing 



by a noise process obtained via sampling at each iterate n a random variable with uniform density 
and intensity a. Experiments with colored, rather than white noise and multiplicative, rather than 
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additive noise led to the same qualitative results. The results for additive white noise are shown in 
Fig. [8] for three different levels of noise intensity a. 

As expected, only the largest tongues survive the presence of the noise; in particular, there is 
only a finite number of surviving tongues, shown in red in Fig. [8j Within such a surviving tongue, 
the random attractor A{uS) is a random periodic cycle of period q (not shown). In the blue region 
outside the Arnol'd tongues, the random attractor is a fixed but random point A{uS) — {a(cj)}: if 
one starts a numerical simulation for a fixed realization of the noise cj, all initial data x converge 
to the same fixed point a, say. 



o=0.05 0=0.10 o=0.15 




Effect of the noise on the PDF of the Arnol'd tongue 1/3 



Figure 8. Arnol'd tongues in the presence of additive noise with different noise 
amplitudes a. Upper panels: Arnol'd tongues for a — 0.05,0.10 and 0.15; lower 
panel: PDF for e = 0.9 and the three a-values in the upper panels: a — 0.05 (red 
curve), a = 0.10 (blue curve), and a = 0.15 (black curve). 

We illustrate this remarkable property in Fig. [9] in the case of a random fixed point, for given 
e and r. The Lyapunov exponent for the three distinct trajectories shown in Fig. [9] is strictly 
negative and the trajectories are exponentially attracted to the single random fixed point a (a;), 
the realization of the driving system 6{uS) being the same for all the trajectories; see Appendices A 
and B. Kaijser [102J provided rigorous results on this type of synchronization phenomenon, but in 
a totally different conceptual setting. Interestingly, as the noise intensity increases, the Lyapunov 
exponent becomes more negative, so that the synchronization occurs even more rapidly, given a 
fixed realization uo. 

This clustering behavior of trajectories with different initial data is in fact well known for flows 
on the circle [103J. In our example, this phenomenon in phase space is related to a smoothing of the 
Devil's staircase in parameter space, the latter cannot be solely explained by the former. Indeed, 
we show in Appendix B that for different irrational numbers and a sufficiently high noise level, the 
corresponding stochastic dynamics are stochastically equivalent, an equivalence that results in the 
smoothing of certain steps of the Devil's staircase. 

As shown in the lower panel of Fig. |8j there is also a direct relationship between the random 
dynamics and the support of the PDF on the circle. For a given noise level, this support can 
either be the union of a finite number of disjoint intervals (red and blue curves) or it can fill the 
whole circle (black curve). The random attractor is, accordingly, either a random periodic orbit, 
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Time 



Figure 9. Synchronization by additive noise: three distinct trajectories (in blue, 
red and black) of = F r ^. u (x n ), with F T ^ given by (B.l); the three trajectories 



start from three initial points on the circle, but are driven by the same realization 
uj of the noise, and thus converge to the same random fixed point a (a;), which is 
moving with time. The parameters are e = 0.5, r = 0.283 and a — 0.3, and the 
corresponding Lyapunov exponent is A ~ —0.0104. 



with the disjoint intervals being visited in succession, or a random fixed point; this PDF behavior 
characterizes the level of the noise needed to destroy a given tongue. 

An exact definition of random fixed point and random periodic orbit is given in Appendix 
B, where we provide a rigorous justification of the numerical results in Figs. [7J [HJ and [9| This 
theoretical analysis helps clarify the interaction between noise and nonlinear dynamics in the context 
of the GCM classification problem we are interested in. 

4. Concluding remarks 

We recall that Section 2 dealt with the natural, interannual and interdecadal variability of 
the ocean's wind-driven circulation. The oceans' internal variability is an important source of 
uncertainty in past-climate reconstructions and future-climate projections \T0\ [TTJ [12] . In Section 
3 and Appendices A and B, we dealt more generally with the problem of structural instability 
as a possible cause for the stubborn tendency of the range of uncertainties in climate change 
projections to increase, rather than diminish over the last three decades [TJ El [3j S] ; see again Fig. 
1. We summarize here the main results of the two sections in succession, and outline several open 
problems. 

The wind-driven double- gyre circulation dominates the near-surface flow in the oceans' mid- 
latitude basins. Particular attention was paid to the North Atlantic and North Pacific, traversed 
by the best-known oceanic jets, namely the Gulf Stream and the Kuroshio Extension (see Fig. 2). 
The wind-driven circulation exhibits very rich internal dynamics and multiscale behavior associated 
with turbulent mesoscales (see Fig. 3) . Aside from the intrinsic interest of this problem in physical 
oceanography, these major oceanic currents help regulate the climate of the adjacent continents, 
while their low- frequency variability affects past, present and future global climate. 

Thanks in part to the systematic use of dynamical systems theory, a comprehensive under- 
standing of simple, barotropic, quasi-geostrophic (QG) models of the double-gyre circulation has 
been achieved over the last two decades, and was reviewed in Section 2 here. In particular, the 
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importance of symmetry-breaking and homoclinic bifurcations (see Fig. 4) in explaining the ob- 
served low-frequency variability has been validated across a wide hierarchy of models, including 
models with much more comprehensive physical formulation, more realistic geometry, and greater 
resolution in the horizontal and vertical [TH [12] . This successive-bifurcation theory also provides 
a self-consistent explanation for the climatically important 7-year and 14-year oscillations of the 
oceanic circulation and the related atmospheric phenomena in and around the North- Atlantic basin 

[IDIiaiHEHlEHlEa 

The next challenge in physical oceanography is to reconcile the points of view of dynamical 
systems theory and statistical mechanics in describing the interaction between the largest scales 
of motion and geostrophic mesoscale turbulence, which is fully captured in baroclinic QG models. 
We emphasize that the complexity of these models of the double-gyre circulation is intermediate 
between high-end GCMs and simple "toy" models; these models offer, therefore, an ideal laboratory 
to test our ideas. In particular, stochastic parametrizations of the rectification process, absent in 
barotropic QG models, could be studied using some of the concepts and tools from RDS theory 
presented here. Note that the RDS approach has already been used in the context of stochastic 
partial differential equations, in particular for showing the existence of random attractors, as well 
as stable, unstable and inertial manifolds. Thus RDS concepts and tools are not restricted to 
finite-dimensional systems [104} \W5\ 1106] . 

In Section 3, we have addressed the range-of-uncertainty problem for IPCC-class GCM sim- 
ulations (see Fig. 1) by considering them as stochastically perturbed dynamical systems. This 
approach is consonant with recent interest for stochastic parametrizations in the high-end modeling- 
and-simulation community [TH [75l [76l [771 EE1 EH] . Rigorous mathematical results from the dynam- 
ical systems literature suggest that — in the absence of stochastic ingredients — GCMs as well as 
simpler models, found on the lower rungs of the modeling hierarchy |17| . are bound to differ from 
each other in their results. 

This sensitivity follows from the fact that, among deterministic dynamical systems, those that are 
hyperbolic are essentially the only ones that are also structurally stable, at least in the C 1 case [621 
[631 EU [65] . Thus, because hyperbolic systems are not dense in the set of smooth deterministic ones 
[67] , we are led to conclude that the topological, structural-stability approach does not guarantee 
deterministic-model robustness, in spite of its many valuable contributions so far. Related issues 
for GCM modeling were emphasized recently by Mitchell [83], Held [M] and McWilliams [85] . 

We have gone one step further and considered model robustness in the presence of stochastic 
terms; such terms could represent either parametrizations of unresolved processes in GCMs or 
stochastic components of natural or anthropogenic forcing, such as volcanic eruptions or fluctuations 
in greenhouse gas or aerosol emissions. Despite the obvious gap between idealized models and 
high-end simulations, we have brought to bear random dynamical systems (RDS) theory [20] on 
the former. 

In this framework, we have considered a robustness criterion that could replace structural sta- 
bility, through the concept of stochastic conjugacy (see Figs. 10 and 11). We have shown, for a 
stochastically perturbed Arnol'd family of circle maps, that noise can enhance model robustness. 
More precisely, this circle map family exhibits structurally stable, as well as structurally unstable 
behavior. When noise is added, the entire family exhibits stochastic structural stability, based on 
the stochastic-conjugacy concept, even in those regions of parameter space where deterministic 
structural instability occurs for vanishing noise (see Figs. [7] and [8]). 

Clearly the hope that noise can smooth the very highly structured pattern of distinct behavior 
types for climate models, across the full hierarchy, has to be tempered by a number of caveats. 
First, serious questions remain at the fundamental, mathematical level about the behavior of non- 
hyperbolic chaotic attractors in the presence of noise [86l EH [88] . Likewise, the case of driving by 
nonergodic noise is being actively studied [107} I108[ 1109] . 
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Second, the presence of certain manifestations of a Devil's staircase has been documented across 
the full hierarchy of ENSO models [TT^ l93| \M 1951 1961 1971 198} 199] . as well as in certain observations 
jTTl [99] . Interestingly, both GCMs and observations only exhibit a few, broad steps of the staircase, 
such as 4 : 1 = 4 yr, 4 : 2 = 2 yr, and 4 : 3 = 16 months. Does this result actually support the 
idea that nature and its detailed models always provide sufficient noise to achieve considerable 
smoothing of the much finer structure apparent in simpler models? Be that as it may, we need 
a much better understanding of how different types of noise — additive and multiplicative, white 
and colored — act across even a partial hierarchy of models, say from the simplest ones, like those 
studied in Section 3, to the intermediate ones considered in Section 2. 

Third, one needs to connect more closely the nature of a stochastic parametrization and its 
effects on the model's behavior in phase-parameter space. As shown in Appendix B, not all types 
of noise are equal with respect to these effects. We are thus left with a rich, and hopefully fruitful, 
set of questions, which we expect to pursue in future work. 
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Appendix A. RDSs and random attractors 

We present here briefly the mathematical concepts and tools of random dynamical systems, 
random attractors and stochastic equivalence. We shall use the concept of pullback attractor 
introduced in Section 3.2.1 to define the closely related notion of a random attractor, but need first 
to define an RDS. We denote by T the set Z, for maps, or R, for flows. Let (X, B) be a measurable 
phase space, and (fi, J 7 , P, (9(t)) t eT) be a metric dynamical system i.e. a flow in the probability 
space (fi, J 7 , P), such that (t,uo) \-> 9{t)uo is measurable and 9{t) : —¥ is measure preserving, 
i.e.,0(t)P = P. 

Let (/9:TxfixX^X, (£, u,x) i— p(t,uo)x, be a mapping with the two following properties: 
(Ri): (p(0,(jj) = Idx, and 

(R2) (the cocycle property): For all s, t E T and all uo E ^, 

<p(t + 5, uo) — (/?(£, 9(s)uo) o <p($, uo). 

If p is measurable, it is called a measurable RDS over 9. If, in addition, X is a topological space 
(respectively a Banach space), and <p satisfies (t,u) \-> (p(t,uo)x continuous (resp. C k , 1 < k < 00) 
for all (£, uo) E T x fi, then p is called a continuous (resp. C k ) RDS over the flow 9. If so, then 

(A.l) (u,x) h+ G(t)(x,uo) := (0(t)u,<p(t,u)x), 

is a (measurable) flow on x X, and is called the skew-product of 9 and (p. In the sequel, we shall 
use the terms "RDS" or "cocycle" synonymously. 

The choice of the so-called driving system 9 is a crucial step in this set-up; it is mostly dictated 
by the fact that the coupling between the stationary driving and the deterministic dynamics should 



respect the time invariance of the former, as illustrated in Fig. 10 The driving system 9 also plays 
an important role in establishing stochastic conjugacy |110j and hence the kind of classification we 
aim at. 

The concept of random attractor is a natural and straightforward extension of the definition of 



pullback attractor (3.2), in which Sell's [81] process is replaced by a cocycle, cf. Fig. 10, and the 
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attractor A now depends on the realization uj of the noise, so that we have a family of random 
attractors A(oj), cf. Fig. 
the noise back to t — » - 



11 



Roughly speaking, for a fixed realization of the noise, one "rewinds" 
oc and lets the experiment evolve (forward in time) towards a possibly 
attracting set A{uS)\ the driving system 9 enables one to do this rewinding without changing the 



statistics, cf. Figs. [10|and 11 




Figure 10. Random dynamical systems (RDS) viewed as a flow on the bundle 
X x ft = "dynamical space" x "probability space." For a given state x and realization 
cj, the RDS (f is such that Q(t)(x,u) = (0(t)u, (p(t,u)x) is a flow on the bundle. 



Fullback attraction toA((x>) 




Figure 11. Schematic diagram of a random attractor A{uS), where uj E ft is a fixed 
realization of the noise. To be attracting, for every set B of X in a family 93 of 
such sets, one must have lim^+oo dist(B(#(— t)u), A{uj)) = with B{6{—t)uS) := 
(p(t,0(—t)(jj)B; to be invariant, one must have (p(t,u)A(uj) = A(0{t)uj). This defini- 
tion depends strongly on 53; see [112] for more details. 



Other notions of attractor can be defined in the stochastic context, in particular based on the 
original SDE; see |lllj or |112j for a discussion on this topic. The present definition, though, will 
serve us well. 

Having defined RDSs and random attractors, we now introduce the notion of stochastic equiv- 
alence or conjugacy, in order to rigourously compare two RDSs; it is defined as follows: two 
cocycles ipi(cj,t) and ^{uo^t) are conjugated if and only if there exists a random homeomorphism 
h E Homeo(X) and an invariant set such that h(u)(0) = and 

(A.2) Lpi{u, t) = /i(6>(t)cj) _1 o (p 2 (u, t) o h{uo). 
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Stochastic equivalence extends classic topological conjugacy to the bundle space Ixft, stating 
that there exists a one-to-one, stochastic change of variables that continuously transforms the phase 
portrait of one sample system in X into that of any other such system. 



Appendix B. Coarse-graining of the circle map family 



We provide here a rigorous justification of the numerical results obtained in Section [3. 2. 2| on the 
topological classification of the family of Arnol'd circle maps in the presence of noise. Consider the 
following random family of diffeomorphisms: 



(B.l) F T ^ ;uj (x) := x + r + cru — esin(27nc) mod 1, 

for x E S 1 , e a real parameter in (0,1), and uj a random parameter distributed in the compact 
interval / = [—1/2,1/2] with fixed distribution v and noise intensity a. We denote by F r ^ e the 
corresponding deterministic family of diffeomorphisms when the noise is switched off, a = 0. 

In the RDS framework, we need to specify the metric dynamical system modeling the noise. We 
choose here the interval al as the base for the probability space and define the flow 9 simply as 
mapping the point uj into its successor in a sequence of realizations of the noise. One could also 
use an irrational rotation on ft for instance; in either case, ergodicity is ensured. 

For the sake of simplicity, we omit for the moment the dependence on r and e. In discrete time, 
with T = Z, we define a map (j) : T x ft x S 1 — >> S 1 , (n, Co>, x) i— »* 0(n, uj)x, such that 



Fen-iu o • • • o Ffj, n > 1, 
Ids 1 j n = 0, 



(B.2) </>(n,u) := 

One can prove easily that this (/) satisfies the cocycle property and is in fact a C°° RDS on S 1 over 



The pair of mappings O := (0,0) is the corresponding skew-product (A.l), and it defines a flow 
on x S 1 by the relation: 



(B.3) (u,x) h+ Q(n)(u,x) := ((9 n cj, <f>(n, u)x). 

A stationary measure m on S 1 under the random diffeomorphism F T ^ ;UJ yields a O-invariant measure 
H := m x z/, i.e. O n /x = /x; explicitly, 

f QxSl f(u,x)n(du,dx) = 
/ nxS i f{9 n u,<j>(n,u)x)ii{dw,dx) 

for all n G T and / G L 1 ^ x S 1 ,//). 

Let us recall the following important proposition |113j concerning the stationary measures ob- 
tained from the random family {F T ^ ;u }. 

Theorem B.l. The random circle diffeomorphism F T ^ ]0J has a unique stationary measure m T:6 . 
The support of m T ^ consists either of q mutually disjoint intervals or of the entire circle S 1 . The 
density function (j) T e is in C 00 ^ 1 ) and depends C°° on r. The invariant measure \i is ergodic. If 
the support of m is connected, then it is mixing and so is fi. 

Mixing for m means that, for any bounded function / : S 1 — > R, and for an arbitrary initial 
point xo E S 1 , E(/(0(n, (jj)xq) tends to J gl f{x)m{dx) as n — > +oc; see |114j for more on random 
attractors and mixing. 

For deterministic diffeomorphisms of the circle, the rotation number measures the average rota- 
tion per iterate of F r e . In the presence of noise, one can still define a rotation number for F r ^ ;uj , 
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namely 

F k . (x) 

(B.5) Pt^uj{x) = lim /e,UJ 



k^oo k 

where F denotes the lift of a map F, acting on S 1 modulo 1, to a map acting on R. For fixed r 
and e, we can then show that p r ^-uj exists for ^/-almost all uj and is a constant; this constant p T]UJ 
is independent of x and uj |113j . Furthermore, r — >► p r ^ e is C°° for each e, which is not true in the 
deterministic case with a = 0; see again |113j . 



Theorem |B.1| has a natural geometric counterpart in terms of random attractors, as confirmed 
through our numerical study; see again Fig. [8j More precisely, we introduce also the following 
definitions of random fixed point and random periodic orbit; these definitions differ somewhat from 
those given in [113J. 

Definition B.l. A random fixed point is a measurable map a : — » S 1 for which 
(B.6) <f>(l,u)a(uj) = a(0(w)), 

/or v -almost all uj E fi, i.e. s^c/i t/&at f2 x a(fi) an invariant set for the flow given by the skew- 
product 0. A random periodic orbit of period q is likewise an invariant set with cardinality q in 
fibers S 1 x {uj} for is -almost all uj. 

With these definitions, the following results of |113j still hold. 

Theorem B.2. For a random diffeomorphism F T ^ U of the circle S 1 , with a stationary measure m 
supported on a union E of q disjoint intervals, the corresponding skew-product restricted to E 
has precisely one attracting random periodic orbit and one repelling random periodic orbit. 

Attraction in the preceding theorem means that lim^oo \F™{x) — F™(cl(uj))\ = 0, for a set of 
initial data (x,uj) E S 1 x Q with positive A x ^/-measure, in the case of a random attracting fixed 
point; here A is Lebesgue measure on S 1 and the extension to a random periodic orbit is obvious. 

Using these two theorems and rigorous results on random point attractors |112| . we can show 
that (i) if the support of the stationary measure is the whole circle (black curve in Fig. [8]), then 
there exists one random fixed point which is pullback attracting; and (ii) if the support consists of 
q disjoint intervals, then the random attractor is a random periodic orbit of period q (red and blue 
curves) . 

Having explained how the connectedness of the PDF support at different noise levels is related 
to the nature of the random attractor, we now turn to an explanation of the "disappearance" 
of the smaller steps in the Devil's staircase, as the noise level increases. To do so, we consider 
the Lyapunov spectrum of an RDS, which still relies on the Oseledets |115j multiplicative ergodic 
theorem (MET). 

To state an MET for RDS on manifolds, we differentiate 0(n,c<j) at x E S 1 , and obtain the linear 
map 

(B.7) T<P(n, uj, x) : T X M -> T^ )x M, 

where T<j) is a continuous linear cocycle on the tangent bundle TM of the manifold M over the 
skew-product flow 0. If the flow (j) possesses an ergodic invariant measure fi such that the required 
integrability condition for applying the MET is verified with respect to /i, then the MET holds for 
(j) over M [TTH] . 



Because of Theorem B.l here, we can apply the MET to our problem and conclude that a unique 
Lyapunov exponent exists for the linearization of each diffeomorphism belonging to our family of 
random diffeomorphisms, and that this exponent is independent of the realization of the noise. We 
show next how to use the Lyapunov spectrum in studying the stochastic equivalence classes of a 
given RDS family, along with its driving system 9. This last aspect of the classification problem is 
outlined for linear hyperbolic cocycles. 
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N.D. Cong |117j has shown that, even in the linear context, the main difference with respect to 
the deterministic case is that the classification depends strongly on the properties of 0, which is 
directly linked to the system noise and its modeling. For instance, if 9 is an irrational rotation on 
S 1 , one can construct infinitely many classes of hyperbolic cocycles that are not pairwise topologi- 
cally equivalent, by playing essentially with the orientations of the cocyles, i.e. reversing between 
clockwise and anticlockwise rotation on S 1 . As we shall see, such difficulties can be avoided in 
the case of noisy Arnol'd tongues, especially for additive noise. Related issues still form an active 
research area in RDS theory; see |116| for a brief survey. 

A key ingredient for the linear classification is the notion of coboundary, which we recall herewith. 

Definition B.2. A measurable set K C is called a coboundary if there exists a set H E T such 
that K — HA9H, where HA9H denotes the symmetric difference of H and 9H. 

Let A and B be two linear random maps on and denote by deg A(uj) and deg B(uo) the 
degrees of the maps A(uj) and B(uj) with respect to a chosen random orientation. These degrees 
are just the sign of the determinant of the corresponding random matrices, and equal —1 or 1; see 
[HO} I117j for details. Consider the two linear hyperbolic cocycles &a and associated with the 
maps A and £>, and the following subset of 

(B.8) Cab = {u E n\deg A(uo) • deg B(u) = -1}; 

Cab is just the set of all uj E £1 for which the degrees of the two linear maps A{uj) and B(uj) differ. 
The main theorem for the classification of our diffeomorphisms of S 1 follows |117j . 

Theorem B.3. Two one- dimensional linear hyperbolic cocycles $^ and <&b are conjugate if and 
only if the following conditions hold: 

(i) sign \a — sign \b, and 

(ii) the associated set Cab is a coboundary. 

Here Xa and Xb indicate the Lyapunov exponents of <&a and respectively. 

Before applying this result, let us explain heuristically how a Devil's staircase step that corre- 
sponds to a rational rotation number can be "destroyed" by a sufficiently intense noise. Consider 
the period-1 locked state in the deterministic setting. At the beginning of this step, a pair of fixed 
points is created, one stable and the other unstable. As the bifurcation parameter is increased, 
these two points move away from each other, until they are tv radians apart. Increasing the param- 
eter further causes the fixed points to continue moving along, until they finally meet again and are 
annihilated in a saddle-node bifurcation, thus signaling the end of the locking interval. 

When noise is added, we have to distinguish between a "strongly locked" regime, where the 
stable and unstable fixed points are nearly tt radians apart, and a "weakly locked" regime, where 
these two fixed points are close to each other. In both regimes, the relaxation time in the vicinity of 
the stable point represents an important time scale of the problem. In the strongly locked regime, 
this is the only time scale of interest. In the weakly locked regime, though, the process of escaping 
across the unstable fixed point is nonnegligible and the associated escape time becomes the second 
time scale of interest. From these heuristic considerations it follows that the distinction between 
strong and weak locking depends on the strength of the external noise. 

If we consider period-T locked states, with T > 2, the same kind of reasoning can be applied to 
the stable and unstable T-cycle. We conclude therefore, for a fixed e > 0, that the narrower Devil's 
staircase steps are the least robust, while the wider ones are the most robust. 

The fact that a locked case becomes unlocked when noise is growing implies in particular that the 
rotation number p r?€ becomes irrational for a sufficiently high noise level. According to Theorem 
2.1 of [102] . the Lyapunov exponent is strictly negative in this case almost surely. Moreover, by 
reinterpreting other results of Kaijser [102J in our RDS framework, we can show that the random 
attractor is in fact a random fixed point; this, in turn, allows us to conclude that the corresponding 
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linearized cocycle at this random fixed point is hyperbolic. Next, by using the Hartman-Grobman 
theorem for RDSs |118( 1119] 1120] . we can conjugate the nonlinear cocycle with its linearization; in 
fact, Theorem 3.1 of |119j says that this conjugacy is global. 

Consider now two linearized cocycles $^ and <E>#, at one and the same or at two distinct random 
fixed points of the family of random diffeomorphisms, for the same noise intensity, and denote by 
A((jj) and B{uj) the random linear parts of the cocycles ^a and <&b respectively; it follows from 
our model of noisy circle maps that Cab is empty. Indeed, the noise being additive, the random 
orientation is preserved for different parameter values. But 9 is assumed to be ergodic, and so we 
have that f2A0f2 is empty and, therefore, Cab is a coboundary. Therewith, Theorem B.3 can be 
applied to obtain the desired result for the problem considered here: with an appropriate amount 
of noise, two deterministic diffeomorphisms that are not topologically equivalent can fall into the 
same topological stochastic class! The numerical results of Section 3.2 are entirely in agreement 
with this assertion. 

Note that the set Cab could differ from a coboundary, if the noise occurred additively in the 
phase of the nonlinear term, for instance. Here we see the importance of noise modeling in obtaining 
the conjectural view of Fig. [6] for a family of dynamical systems in general. 

It follows, in particular, that the exact nature of the stochastic parametrizations in a family 
of GCMs does matter. It's not enough to follow the trend by devising and implementing such 
parametrizations: one should test that a given parametrization, once found to be suitable in other 
respects, does improve the proximity, in an appropriate sense, between climate simulations within 
the family of GCMs for which it has been been developed. 
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